---
  title: "CARRS RDD analysis :Robustness check keeping only follow uo observations within 12 months"
knit: (function(input, ...) {rmarkdown::render(input, output_file="")})
output: pdf_document
geometry: "left=2cm,right=2cm,top=1cm,bottom=1cm"
date: "`r format(Sys.time(), '%d %B, %Y')`"
---

```{r , setup, include=FALSE, eval = TRUE}
knitr::opts_chunk$set(
   comment = '', warning=FALSE,  message=FALSE, results= FALSE , fig.width = 9,fig.height = 7, echo=FALSE, root.dir =""
)

 knitr::opts_knit$set(
   root.dir = ""    
 )
```


```{r packages}
rm(list=ls())
library(tidyverse)
library(gridExtra)
library(dplyr)
library(rdrobust)
library(rddensity)
library(viridisLite)
library(foreign)
library(ggplot2)
library(data.table)
library(doBy)
library(socviz)
library(openxlsx) 
library(modelsummary)
library(labelled)
library(gtsummary)
library(gt)
library(kableExtra)
```

```{r loaddata}
data_treatment<-data.frame(read.csv("carrs_full.csv"))
data_treatment<-as.data.table(data_treatment)

data_diagnosis<-data.frame(read.csv("carrs_selected.csv"))
data_diagnosis<-as.data.table(data_diagnosis)

data_bp<-data.frame(read.csv("carrs_bp.csv"))
data_bp<-as.data.table(data_bp)
```
##calculating the bandwidth for each sub group

```{r}
# Males and females separately
rv.update = function(b,c) {
  out=c[female==b,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                            (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
  a_b<<-out
}

rv.update(1,data_diagnosis)
diag_female<-a_b[female==1,]
rm(a_b)

rv.update(0,data_diagnosis)
diag_male<-a_b[female==0,]
rm(a_b)

rv.update(1,data_treatment)
treat_female<-a_b[female==1,]
rm(a_b)

rv.update(0,data_treatment)
treat_male<-a_b[female==0,]
rm(a_b)

# Total population
rv.updateMT = function(c) {
  out=c[,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                             (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
  a_b<<-out
}

rv.updateMT(data_diagnosis)
diag_total<-a_b
rm(a_b)

rv.updateMT(data_treatment)
treat_total<-a_b
rm(a_b)
```
##RDD Outputs function

```{r fuc}
carrs_rdd = function(a,b,d) {
  d$Y<-a
  d$X<-b
 
  out =  rdrobust(y = d$Y, d$X, kernel = "triangular", p = 2, q = 3, covs = cbind(d$w01, d$w23, d$w45),
                  c = 0, cluster = d$pid, rho = 1 , all = T)
  summary(out)
  
  ret <- data.frame(
    sex = c(as.character(round(100*out$coef[1,1],2)), as.character(round(out$pv[3,1],2)), 
            paste0( "(", as.character(round(100*out$ci[3,1],2)), " - ", as.character(round(100*out$ci[3,2],2)), ")" ), 
            as.character(round(out$bws[1,1],2)), as.character(out$N_h[1] + out$N_h[2]))
  )
  rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  
  output<<-ret %>% add_rownames  
}
```

##Calling the specific subroups and saving the output

```{r diagnosis}
rdd = function(a) {
  carrs_rdd(a$hypdiag,a$runningvar,a)
}

rdd(diag_total)
diag_total <-output%>%mutate(Total=sex)%>%select(-sex)

rdd(diag_male)
diag_male <-output%>%mutate(Male=sex)%>%select(-sex)

rdd(diag_female)
diag_female <-output%>%mutate(Female=sex)%>%select(-sex)
```

```{r treatment}
rdd = function(a) {
  carrs_rdd(a$hypdrug,a$runningvar,a)
}

rdd(treat_total)
treat_total <-output%>%mutate(Total=sex)%>%select(-sex)

rdd(treat_male)
treat_male <-output%>%mutate(Male=sex)%>%select(-sex)

rdd(treat_female)
treat_female <-output%>%mutate(Female=sex)%>%select(-sex)
```

#combining males and females columsn

```{r}
total_table<-merge(diag_total,treat_total, by="rowname")
total_table<-total_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))


colnames(total_table)[2] <- "Diagnosis"   
colnames(total_table)[3] <- "Treatment"   

female_table<-merge(diag_female,treat_female, by="rowname")
female_table<-female_table%>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(female_table)[2] <- "Diagnosis"   
colnames(female_table)[3] <- "Treatment"    

male_table<-merge(diag_male,treat_male, by="rowname")
male_table<-male_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(male_table)[2] <- "Diagnosis" 
colnames(male_table)[3] <- "Treatment"    
```

##table making function
```{r }

  filename <- paste0("quad_combined.tex")
  caption=paste0("Effect of home-based screening on hypertension diagnosis and treatment: Local quadratic trend")

comb_table<-rbind(total_table,female_table,male_table)

combined<-kbl(comb_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lcccc", format = "latex",col.names = c(" ","Diagnosis","Treatment"),digits = 2) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10) %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15) %>%
    footnote(general =c("The model estimation included a local quadratic trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and CIs resulting from robust bias-corrected standard errors clustered at the individual-level.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")
```

##### Difference in BP level #####


##calculating the bandwidth for each sub group

```{r RVdiff}
rv.updateBP = function(b,c) {
  out=c[female==b,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                            (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
  a_b<<-out
}

rv.updateBP(1,data_bp)
BP_female<-a_b[female==1,]
rm(a_b)

rv.updateBP(0,data_bp)
BP_male<-a_b[female==0,]
rm(a_b)

# Total population
rv.updateMTBP = function(c) {
  out=c[,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                             (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
  a_b<<-out
}

rv.updateMTBP(data_bp)
BP_total<-a_b
rm(a_b)

```
##RDD Outputs function
```{r fucDiff}
carrs_rddBP = function(a,b,d) {
  d$Y<-a
  d$X<-b
 
  out =  rdrobust(y = d$Y, d$X, kernel = "triangular", p = 2, q = 3, covs = cbind(d$w02, d$w24),
                  c = 0, cluster = d$pid, rho = 1)
  summary(out)
  
  ret <- data.frame(
    sex = c(as.character(round(out$coef[1,1],2)), as.character(round(out$pv[3,1],2)), 
            paste0( "(", as.character(round(out$ci[3,1],2)), " - ", as.character(round(out$ci[3,2],2)), ")" ), 
            as.character(round(out$bws[1,1],2)), as.character(out$N_h[1] + out$N_h[2]))
  )
  rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  
  output<<-ret %>% add_rownames  
}
```

##Calling the specific subroups and saving the output

```{r syssdiff}
rddBP = function(a) {
  carrs_rddBP(a$sysdiff,a$runningvar,a)
}

rddBP(BP_total)
sysdiff_total <-output%>%mutate(Total=sex)%>%select(-sex)

rddBP(BP_male)
sysdiff_male <-output%>%mutate(Male=sex)%>%select(-sex)

rddBP(BP_female)
sysdiff_female <-output%>%mutate(Female=sex)%>%select(-sex)
```

```{r diasdiff}
rddBP = function(a) {
  carrs_rddBP(a$diasdiff,a$runningvar,a)
}

rddBP(BP_total)
diasdiff_total <-output%>%mutate(Total=sex)%>%select(-sex)

rddBP(BP_male)
diasdiff_male <-output%>%mutate(Male=sex)%>%select(-sex)

rddBP(BP_female)
diasdiff_female <-output%>%mutate(Female=sex)%>%select(-sex)
```

#combining males and females columsn

```{r}
totalBP_table<-merge(sysdiff_total,diasdiff_total, by="rowname")
totalBP_table<-totalBP_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))


colnames(totalBP_table)[2] <- "Systolic"   
colnames(totalBP_table)[3] <- "Diastolic"   

femaleBP_table<-merge(sysdiff_female,diasdiff_female, by="rowname")
femaleBP_table<-femaleBP_table%>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(femaleBP_table)[2] <- "Systolic"   
colnames(femaleBP_table)[3] <- "Diastolic"    

maleBP_table<-merge(sysdiff_male,diasdiff_male, by="rowname")
maleBP_table<-maleBP_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(maleBP_table)[2] <- "Systolic" 
colnames(maleBP_table)[3] <- "Diastolic"    
```

##table making function
```{r }

  filename <- paste0("quad_BPdiff_combined.tex")
  caption=paste0("Effect of home-based screening on change in systolic and diastolic blood pressure: Local quadratic trend")

combBP_table<-rbind(totalBP_table,femaleBP_table,maleBP_table)

combined<-kbl(combBP_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lcccc", format = "latex",col.names = c(" ","Systolic","Diastolic"),digits = 2) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10) %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15) %>%
    footnote(general =c("Change in BP"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")

```

